#ifndef WAVLET
#define WAVLET

#include<iostream>
#include<cmath>
#include<stdio.h>
#include<stdlib.h>
#include <vector>

class Wavelet
{
private:

public:
  virtual double** WT()=0;
  virtual double** W()=0;
  virtual void init(int size, bool inverse)=0;
  virtual void reset()=0;
  virtual int H()=0;
};

class LeGall : public Wavelet
{
private:
    double** w;
    double** wt;
    int size;
    int h;

public:
  
  
    ~LeGall()
    {
        reset();
    }
  
    LeGall(int s, bool inverse):size(s)
    {
	h=8;
	init(s,inverse);
    }
  
  
    void reset()
    {
        for(int i=0; i < size; ++i)
        {
            delete [] w[i];
	    delete [] wt[i];
        }
        delete [] w;
	delete [] wt;
    }
    
    void init(int s, bool inverse)
    {
	size=s;
      
        w= new double* [size];
        wt= new double* [size];


        for(int i=0; i<size; ++i)
        {
            w[i]=new double[size];
            wt[i]=new double[size];
        }
	
      
	double hh[size/2][size];
	double gg[size/2][size];
      
	
		
	for(int i=0; i<size/2; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                hh[i][j]=0;
                gg[i][j]=0;
            }
        }
      
	if(!inverse)
	{      
	    double h[5];
	    double g[3];
	     
	    
	    
	/*    
	waveFun[1]=0.5
        waveFun[2]=-1./std::sqrt(2);
        waveFun[3]=0.5/std::sqrt(2);
        
	
        scaleFun[0]=(-1./8.)*std::sqrt(2);
        scaleFun[1]=(2./8.)*std::sqrt(2);
        scaleFun[2]=(6./8.)*std::sqrt(2);
        scaleFun[3]=(2./8.)*std::sqrt(2);
        scaleFun[4]=(-1./8.)*std::sqrt(2);*/
	    
	    
	    
	    h[0] = (-1.0 / 8.0)*std::sqrt(2.0);
	    h[1] = (2.0 / 8.0)*std::sqrt(2.0);
	    h[2] = (6.0 / 8.0)*std::sqrt(2.0);
	    h[3] = (2.0 / 8.0)*std::sqrt(2.0);
	    h[4] = (-1.0 / 8.0)*std::sqrt(2.0);
      
	    g[0] = -0.5/std::sqrt(2.0);//-1.0 / 2.0;
	    g[1] = 1./std::sqrt(2.0);//1.0;
	    g[2] = -0.5/std::sqrt(2.0);//-1.0 / 2.0; 
	    
	    int j=0;
	    for(int i=0; i<size/2; ++i)
	    {
		    for(int k=0; k<5; k++)
		    {
			if(j+k-2>=0&&j+k-2<=size-1)
			{
			  hh[i][j+k-2]=h[k];
			}else if(j+k-2<0)
			{
			  hh[i][size+k-2]=h[k];
			}else if(j+k-2>size-1)
			{
			  hh[i][size-4-j+k-2]=h[k];
			}
		    }
	      
		    for(int k=0; k<3; k++)
		    {
			if(j+k-0>=0&&j+k-0<=size-1)
			{
			  gg[i][j+k-0]=g[k];
			}else if(j+k-0<0)
			{
			  gg[i][size+j+k-0]=g[k];
			}else if(j+k-0>size-1)
			{
			  gg[i][size-4-j+k-0]=g[k];
			}
		    }
		    j+=2;
	    }
	}else
	{
	      double hi[5];
	      double gi[3];
	    	      
	      
	      hi[0] = (1.0 / 2.0)/std::sqrt(2.0);
	      hi[1] = 1./std::sqrt(2.0);
	      hi[2] = (1.0 / 2.0)/std::sqrt(2.0);

	      gi[0] = (-1.0 / 8.0)*std::sqrt(2.0);//-
	      gi[1] = (-2.0 / 8.0)*std::sqrt(2.0);//-
	      gi[2] = (6.0 / 8.0)*std::sqrt(2.0);
	      gi[3] = (-2.0 / 8.0)*std::sqrt(2.0);//-
	      gi[4] = (-1.0 /8.0)*std::sqrt(2.0);//-

	    int j=0;
	    for(int i=0; i<size/2; ++i)
	    {
		    for(int k=0; k<3; k++)
		    {
			if(j+k-1>=0&&j+k-1<=size-1)
			{
			  hh[i][j+k-1]=hi[k];
			}else if(j+k-1<0)
			{
			  hh[i][size+j+k-1]=hi[k];
			}else  if(j+k-1>size-1)
			{
			  hh[i][size-2-j+k-1]=hi[k];
			}
		    }
	      
		    for(int k=0; k<5; k++)
		    {
			if(j+k-1>=0&&j+k-1<=size-1)
			{
			  gg[i][j+k-1]=gi[k];
			}else if(j+k-1<0)
			{
			  gg[i][size+j+k-1]=gi[k];
			}else  if(j+k-1>size-1)
			{
			  gg[i][size-4-j+k-1]=gi[k];
			}
		    }
		    j+=2;
	    }
	
	      
	}
        
        
        for(int i=0; i<size/2; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                w[i][j]=hh[i][j];
                w[i+size/2][j]=gg[i][j];
            }
        }

        for(int i=0; i<size; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                wt[i][j]=w[j][i];
            }
        }
      
      
    }
    
    
  
  
    double** WT()
    {
        return wt;
    }

    double** W()
    {
        return w;
    }
    
    int H()
    {
      return h; 
    }
    
  
};


class Daub97 : public Wavelet
{
private:
    double** w;
    double** wt;
    int size;
    int h;

public:
  
  
    ~Daub97()
    {
        reset();
    }
  
    Daub97(int s, bool inverse):size(s)
    {
	h=6;
	init(s,inverse);
    }
  
  
    void reset()
    {
        for(int i=0; i < size; ++i)
        {
            delete [] w[i];
	    delete [] wt[i];
        }
        delete [] w;
	delete [] wt;
    }
    
    void init(int s, bool inverse)
    {
	size=s;
      
        w= new double* [size];
        wt= new double* [size];


        for(int i=0; i<size; ++i)
        {
            w[i]=new double[size];
            wt[i]=new double[size];
        }
	
      
	double hh[size/2][size];
	double gg[size/2][size];
      
	
		
	for(int i=0; i<size/2; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                hh[i][j]=0;
                gg[i][j]=0;
            }
        }
      
	if(!inverse)
	{      
	    double h[9      ];
	    double g[7];
	     
	    
      h[0] = 0.026729;
      h[1] = -0.016864;
      h[2] = -0.078223;
      h[3] = 0.266864;
      h[4] = 0.602949;
      h[5] = 0.266864;
      h[6] = -0.078223;
      h[7] = -0.016864;
      h[8] = 0.026729;

      //Filter g0(7, 2);
      g[0] = 0.091272;
      g[1] = -0.057544;
      g[2] = -0.591272;
      g[3] = 1.115087;
      g[4] = -0.591272;
      g[5] = -0.057544;
      g[6] = 0.091272;
	

	    
	    int j=0;
	    for(int i=0; i<size/2; ++i)
	    {
		    for(int k=0; k<9; k++)
		    {
			if(j+k-4>=0&&j+k-4<=size-1)
			{
			  hh[i][j+k-4]=h[k];
			}else if(j+k-4<0)
			{
			  hh[i][size+k-4]=h[k];
			}else if(j+k-4>size-1)
			{
			  hh[i][size-4-j+k-4]=h[k];
			}
		    }
	      
		    for(int k=0; k<7; k++)
		    {
			if(j+k-2>=0&&j+k-2<=size-1)
			{
			  gg[i][j+k-2]=g[k];
			}else if(j+k-2<0)
			{
			  gg[i][size+j+k-2]=g[k];
			}else if(j+k-2>size-1)
			{
			  gg[i][size-4-j+k-2]=g[k];
			}
		    }
		    j+=2;
	    }
	}else
	{
	      double hi[5];
	      double gi[3];
	    	      
	      
	      hi[0] = (1.0 / 2.0)/std::sqrt(2.0);
	      hi[1] = 1./std::sqrt(2.0);
	      hi[2] = (1.0 / 2.0)/std::sqrt(2.0);

	      gi[0] = (-1.0 / 8.0)*std::sqrt(2.0);//-
	      gi[1] = (-2.0 / 8.0)*std::sqrt(2.0);//-
	      gi[2] = (6.0 / 8.0)*std::sqrt(2.0);
	      gi[3] = (-2.0 / 8.0)*std::sqrt(2.0);//-
	      gi[4] = (-1.0 /8.0)*std::sqrt(2.0);//-

	    int j=0;
	    for(int i=0; i<size/2; ++i)
	    {
		    for(int k=0; k<3; k++)
		    {
			if(j+k-1>=0&&j+k-1<=size-1)
			{
			  hh[i][j+k-1]=hi[k];
			}else if(j+k-1<0)
			{
			  hh[i][size+j+k-1]=hi[k];
			}else  if(j+k-1>size-1)
			{
			  hh[i][size-2-j+k-1]=hi[k];
			}
		    }
	      
		    for(int k=0; k<5; k++)
		    {
			if(j+k-1>=0&&j+k-1<=size-1)
			{
			  gg[i][j+k-1]=gi[k];
			}else if(j+k-1<0)
			{
			  gg[i][size+j+k-1]=gi[k];
			}else  if(j+k-1>size-1)
			{
			  gg[i][size-4-j+k-1]=gi[k];
			}
		    }
		    j+=2;
	    }
	
	      
	}
        
        
        for(int i=0; i<size/2; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                w[i][j]=hh[i][j];
                w[i+size/2][j]=gg[i][j];
            }
        }

        for(int i=0; i<size; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                wt[i][j]=w[j][i];
            }
        }
      
      
    }
    
    
  
  
    double** WT()
    {
        return wt;
    }

    double** W()
    {
        return w;
    }
    
    int H()
    {
      return h; 
    }
    
  
};





class NewHarr : public Wavelet
{
private:

    double** w;
    double** wt;

    int size;
    int h;

public:
  
  
    
    ~NewHarr()
    {
	reset();
    }
   
   
    void reset()
    {
	for(int i=0; i < size; ++i)
        {
            delete [] w[i];
	    delete [] wt[i];
        }
        delete [] w;
	delete [] wt;
    }
  
    NewHarr(int s):size(s)
    {
        init(s,true);
	h=2;
    }
    

    int H()
    {
      return h; 
    }
    
    
    
    void init(int s, bool inverse)
    {
      size=s;
      
      double c = std::sqrt(2.0);

        double h [size/2][size];
        double g [size/2][size];

        w= new double* [size];
        wt= new double* [size];


        for(int i=0; i<size; ++i)
        {
            w[i]=new double[size];
            wt[i]=new double[size];
        }

        for(int i=0; i<size/2; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                h[i][j]=0;
                g[i][j]=0;
            }
        }

	int j=0;
        for(int i=0; i<size/2; ++i)
        {
                h[i][j]=c/2.0;
                h[i][j+1]=c/2.0;
                g[i][j]=-c/2.0;
                g[i][j+1]=c/2.0;
                j+=2;
        }

        for(int i=0; i<size/2; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                w[i][j]=h[i][j];
                w[i+size/2][j]=g[i][j];
            }
        }
        
        for(int i=0; i<size; ++i)
        {
            for(int j=0; j<size; ++j)
            {
                wt[i][j]=w[j][i];
            }
        }
      
    }
   
    double** WT()
    {
	return wt;
    }

    double** W()
    {
        return w;
    }
};

class WaveletTransform
{
private:
    Wavelet *wavelet;
public:
    WaveletTransform(Wavelet *w):wavelet(w){}
    ~WaveletTransform()
    {
        delete wavelet;
    }

    void multiplication(double **A, double **B,double**C, int size);
    double** transform(double **Tab,int width, int height);
    double** reverseTransform(double **transformTab,int width, int height);
};


/*
class Harr: public Wavelet
{
private:
    double scaleFun[2];
    double waveFun[2];
    double waveLength;
public:
    Harr()
    {
        waveFun[0]=1./1.4142135623730951;
        waveFun[1]=-1./1.4142135623730951;
        scaleFun[0]=1./1.4142135623730951;
        scaleFun[1]=1./1.4142135623730951;
        waveLength=2;
    }
    double* ScalFunction()
    {
        return scaleFun;
    }
    double* WavelFunction()
    {
        return waveFun;
    }

    double getWaveLength()
    {
        return waveLength;
    }


    double* transform(double *tab, int length)
    {
        double *transformTab=new double[length];

        for(int i=0; i<length; ++i)
        {
            transformTab[i]=0;
        }

        int h = length >> 1;

        int k=0;

        for(int i=0; i<h; ++i)
        {
            for( int j = 0; j < waveLength; ++j )
            {
                k = ( i << 1 ) + j;
                while( k >= length )
                    k -= length;
                transformTab[ i ] += tab[ k ] * scaleFun[ j ];
                transformTab[ i + h ] +=tab[ k ] * waveFun[ j ];
            }
        }
        return transformTab;
    }



    double* reverse( double* transformTab, int length )
    {
        double *reverseTab=new double[length];
	
        for(int i=0; i<length; ++i)
        {
            reverseTab[i]=0;
        }

        int h = length >> 1;

        int k=0;

        for( int i = 0; i < h; i++ )
        {
            for( int j = 0; j < waveLength; j++ )
            {
                k = ( i << 1 ) + j;
                while( k >= length )
                    k -= length;
                reverseTab[ k ] += transformTab[ i ] * scaleFun[ j ] + transformTab[ i + h ] * waveFun[ j ];
                //reverseTab[ k ] *= 0.5;
            }
        }
        return reverseTab;
    }
};
*/


#endif
